%% AFI B1 Mapping Script
%

%% Prepare Workspace
% Clear all variables, and clear the command window.

clear
clc

%% Set the AFI input sequence parameters.
% Note that TR1 = N*TR2.

N=5;                % TRratio
NominalAngle=60;    % Nominal flip angle (degrees)

%% Load the minc images into the Workspace
% Note that using the "loadminc" inputs the data into an in32 format.

afi1 = loadminc('boudreau_vfa_20120314_20120314_094646_24d1_mri.mnc');
afi2 = loadminc('boudreau_vfa_20120314_20120314_094646_25d2_mri.mnc');

afi3 = loadminc('boudreau_vfa_20120314_20120314_094646_26d1_mri.mnc');
afi4 = loadminc('boudreau_vfa_20120314_20120314_094646_27d2_mri.mnc');

%% Call AFIB1Map function.
% This function calculates the B1 maps using two AFI images, TRratio (N)
% and the nominal flip angle (NominalAngle)

b1map=AFIB1Map(afi1,afi2,N,NominalAngle);

b1maphires=AFIB1Map(afi3,afi4,N,NominalAngle);

%% B1 map plots
% 

% Plot histogram
reshapedb1map= reshape(b1map,1, length(b1map(:,1,1))*length(b1map(1,:,1))*length(b1map(1,1,:)));
reshapedb1map(reshapedb1map==0|reshapedb1map>1.5)=[];
figure()
[y x]=hist(reshapedb1map,100);
normTerm = length(reshapedb1map(:,1,1))*length(reshapedb1map(1,:,1))*length(reshapedb1map(1,1,:))
plot(x,y./normTerm)

hold on

reshapedb1maphires= reshape(b1maphires,1, length(b1maphires(:,1,1))*length(b1maphires(1,:,1))*length(b1maphires(1,1,:)));
reshapedb1maphires(reshapedb1maphires==0|reshapedb1maphires>1.5)=[];
[y x]=hist(reshapedb1maphires,100);
normTerm = length(reshapedb1maphires(:,1,1))*length(reshapedb1maphires(1,:,1))*length(reshapedb1maphires(1,1,:));
plot(x,y./normTerm, 'r')

title('Whole-volume B1 map histogram from AFI comparing a 2mm x 2mm x 5mm and 1mm x 1mm x 1mm acquisitions')
xlabel('B1 multiplicative correction factor (a.u.)')
legend('2mm x 2mm x 5mm', '1mm x 1mm x 1mm')

%% B1 map (2x2x5) plots
% 

figure()
imagesc(b1map(:,:,length(b1map(1,1,:))/2));
title('3T AFI B1 map - 2mm x 2mm x 5mm - zCenter')
colorbar
xlabel('xDirection')
ylabel('yDirection')

figure()
imagesc(squeeze(b1map(:,length(b1map(1,:,1))/2,:)));
title('3T AFI B1 map - 2mm x 2mm x 5mm - yCenter')
colorbar
xlabel('xDirection')
ylabel('zDirection')


figure()
imagesc(squeeze(b1map(length(b1map(:,1,1))/2,:,:)));
title('3T AFI B1 map - 2mm x 2mm x 5mm - xCenter')
colorbar
xlabel('yDirection')
ylabel('zDirection')

%% B1 map (1x1x1) plots
% 

figure()
imagesc(b1maphires(:,:,length(b1maphires(1,1,:))/2));
title('3T AFI B1 map - 1mm x 1mm x 1mm - zCenter')
colorbar
xlabel('xDirection')
ylabel('yDirection')

figure()
imagesc(squeeze(b1maphires(:,length(b1maphires(1,:,1))/2,:)));
title('3T AFI B1 map - 1mm x 1mm x 1mm - yCenter')
colorbar
xlabel('xDirection')
ylabel('zDirection')


figure()
imagesc(squeeze(b1maphires(length(b1maphires(:,1,1))/2,:,:)));
title('3T AFI B1 map - 1mm x 1mm x 1mm - xCenter')
colorbar
xlabel('yDirection')
ylabel('zDirection')